///////////////////////////////////////////////////////////////////////////////////////////////////
/// \file pcps.h
///
/// \brief   the primary header file for the PCPS program.
///
///   In all cpp files, this header file should be the first header file included
///   to ensure that the defninition of NDEBUG precedes all other includes.
///
///////////////////////////////////////////////////////////////////////////////////////////////////

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \mainpage PCPS Program Documentation
///
///////////////////////////////////////////////////////////////////////////////////////////////////

#ifndef PROJECTED_CPS_HEADER
#define PROJECTED_CPS_HEADER

#define NDEBUG

#include <iostream>
#include <fstream>
#include <ctime>
#include <cmath>
#include <set>
#include <map>
#include <list>
#include <vector>
#include <string>
#include <typeinfo>
#include <cassert>
#include <complex>
#include <cstdlib>
#include <cstring>
#include <algorithm>
#include <numeric>
#include <sstream>

#include <mpi.h>
#include <pthread.h>

#include <boost/format.hpp>
#include <boost/shared_ptr.hpp>
#include <boost/shared_array.hpp>

#include <src/exception.h>

///////////////////////////////////////////////////////////////////////////////////////////////////
/// \brief   namespace for most functions and classes in the PCPS program
///
///////////////////////////////////////////////////////////////////////////////////////////////////

namespace pcps {

  // incomplete declarations of some classes
  class Input;
  class Cluster;
  class Arena;
  template <class SCALAR> class OpBase;
  template <class SCALAR> class Hamiltonian;
  template <class SCALAR> class Correlator;
  template <class SCALAR> class Wavefunction;

  /// \brief   empty class used for selecting templates for fermions
  class Part_Fermion {};

  /// \brief   empty class used for selecting templates for spinless fermions
  class Part_SpinlessFermion {};

  /// \brief   empty class used for selecting templates for spins
  class Part_Spin {};

  /// \brief   empty class used for selecting templates for the uniform reference
  class Ref_Uniform {};

  /// \brief   empty class used for selecting templates for the RHF reference
  class Ref_RHF {};

  /// \brief   empty class used for selecting templates for the AGP reference
  class Ref_AGP {};

  /// \brief   empty class used for selecting templates for the pairing reference
  class Ref_Pairing {};

  /// \brief   empty class used for selecting templates for the generalized pairing reference
  class Ref_Pfaffian {};

  //-------------------------------------------------------------------------------
  // pcps::init_mutexes -- initializes the global mutexes
  //-------------------------------------------------------------------------------

  void init_mutexes();

  //-------------------------------------------------------------------------------
  // pcps::destroy_mutexes -- destroys the global mutexes
  //-------------------------------------------------------------------------------

  void destroy_mutexes();

  //-------------------------------------------------------------------------------
  // pcps::alloc_mutex -- returns a pointer to a mutex used by the Arena allocator
  //                      to ensure that its allocation/deallocation is thread safe
  //-------------------------------------------------------------------------------

  pthread_mutex_t * alloc_mutex();

  //-------------------------------------------------------------------------------
  // pcps::rand_mutex -- returns a pointer to a mutex for ensuring that threads
  //                     can safely use the global random number generator
  //-------------------------------------------------------------------------------

  pthread_mutex_t * rand_mutex();

  //-------------------------------------------------------------------------------
  // pcps::Stopwatch -- A class for keeping track of elapsed time
  //-------------------------------------------------------------------------------

  class Stopwatch {

    // constructors
    public:
      Stopwatch() : _running(false), _start_time(0.0), _elapsed_time(0.0) {}

    // public member functions
    public:
      void reset() {
        if (_running)
          throw pcps::Exception("cannot reset Stopwatch because it is running");
        _elapsed_time = 0.0;
      }
      void start() {
        if (_running)
          throw pcps::Exception("cannot start Stopwatch because it is already running");
        _start_time = MPI::Wtime();
        _running = true;
      }
      void stop() {
        if (!_running)
          throw pcps::Exception("cannot stop Stopwatch because it is not running");
        _elapsed_time = _elapsed_time + MPI::Wtime() - _start_time;
        _running = false;
      }
      double elapsed_seconds() const {
        if (_running)
          return (_elapsed_time + MPI::Wtime() - _start_time);
        return _elapsed_time;
      }

    // private data
    private:
      bool _running;
      double _start_time; // in seconds
      double _elapsed_time; // in seconds

    // disabled functions
    private:
      Stopwatch(const Stopwatch &); // disable copy constructor
      Stopwatch &operator=(const Stopwatch &); // disable assignment operator
  };

  //--------------------------------------------------------------------------------
  // pcps::split_string -- Splits a string according to the supplied deliminators
  //--------------------------------------------------------------------------------

  std::vector<std::string> split_string(const std::string &str, const std::string &deliminators = " ");

  //--------------------------------------------------------------------------------
  // pcps::pow -- computes integer powers.  Assumes b >= 0 and returns a^b
  //--------------------------------------------------------------------------------

  template <class T> inline T pow(const T a, int b) {
    assert(b >= 0);
    T retval = 1;
    for ( ; b > 0; b--) retval *= a;
    return retval;
  };

  //--------------------------------------------------------------------------------
  // pcps::n_choose_m -- computes the binomial coefficient
  //--------------------------------------------------------------------------------

  inline double n_choose_m(int n, int m) {
    if (n < 0 || m < 0 || m > n) return 0.0;
    double retval = 1.0;
    while (m > 0) retval = ( retval * (n--) ) / (m--);
    return retval;
  }

  //-------------------------------------------------------------------------------
  // pcps::n_integer_solutions -- Computes the number of solutions to the equation
  //                              x(1) + x(2) + ... + x(n) = r
  //                              when the variables x(i) are constrained to be
  //                              integers in the range (0, 1, 2, ..., k)
  //-------------------------------------------------------------------------------

  double n_integer_solutions(const int n, const int r, const int k);

  //-------------------------------------------------------------------------------
  // pcps::get_all_site_pairs -- returns a set containing all pairs of sites
  //-------------------------------------------------------------------------------

  std::set< std::pair<int,int> > get_all_site_pairs(const int nsites);

  //-------------------------------------------------------------------------------
  // pcps::get_nn_site_pairs -- returns a set containing all nearest neighbor pairs
  //                            of sites
  //-------------------------------------------------------------------------------

  std::set< std::pair<int,int> > get_nn_site_pairs(const pcps::Input &userinp);

  //--------------------------------------------------------------------------------
  // pcps::driver -- controls which calculations to perform
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void driver(const pcps::Input & userinp);

  //--------------------------------------------------------------------------------
  // pcps::bfgs_driver -- performs a variational bfgs optimization of the
  //                      wavefunction
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void bfgs_driver(const pcps::Input & userinp,
                   const pcps::Hamiltonian<SCALAR> & ham,
                   pcps::Wavefunction<SCALAR> & wfn);

  //--------------------------------------------------------------------------------
  // pcps::projected_driver -- performs a projected cps optimization of the
  //                           wavefunction
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void projected_driver(const pcps::Input & userinp,
                        const pcps::Hamiltonian<SCALAR> & ham,
                        pcps::Wavefunction<SCALAR> & wfn);

  //--------------------------------------------------------------------------------
  // pcps::imag_time_driver -- performs an imaginary time optimization of the
  //                           wavefunction
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void imag_time_driver(const pcps::Input & userinp,
                        const pcps::Hamiltonian<SCALAR> & ham,
                        pcps::Wavefunction<SCALAR> & wfn);

  //--------------------------------------------------------------------------------
  // pcps::energy -- computes the energy of the cps wavefunction
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  SCALAR energy(const pcps::Input & userinp,
                const pcps::Hamiltonian<SCALAR> & ham,
                const pcps::Wavefunction<SCALAR> & wfn,
                const int perturbation_level);

  //-------------------------------------------------------------------------------
  // pcps::steepest_descent -- updates the correlators using a steepest descent
  //                           step, where the gradient is evaluated by monte carlo
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void steepest_descent(const pcps::Input & userinp,
                        const pcps::Hamiltonian<SCALAR> & ham,
                        pcps::Wavefunction<SCALAR> & wfn);

  //-------------------------------------------------------------------------------
  // pcps::steepest_descent_detail -- updates the correlators using a steepest
  //                                  descent step, where the gradient is evaluated
  //                                  by monte carlo
  //-------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF>
  void steepest_descent_detail(const pcps::Input & userinp,
                               const pcps::Hamiltonian<SCALAR> & ham,
                               pcps::Wavefunction<SCALAR> & wfn);

  //-------------------------------------------------------------------------------
  // pcps::stoch_imag_time -- Updates the correlators using an imaginary time
  //                          evolution in which each update step is based on
  //                          stochastic information on the wavefunction and its
  //                          derivatives.
  //                          This method is very similar to the stochastic
  //                          reconfiguration method.
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void stoch_imag_time(const pcps::Input & userinp,
                       const pcps::Hamiltonian<SCALAR> & ham,
                       pcps::Wavefunction<SCALAR> & wfn);

  //-------------------------------------------------------------------------------
  // pcps::stoch_imag_time_detail -- Updates the correlators using an imaginary
  //                                 time evolution in which each update step is
  //                                 based on stochastic information on the
  //                                 wavefunction and its derivatives.
  //                                 This method is very similar to the stochastic
  //                                 reconfiguration method.
  //-------------------------------------------------------------------------------

  template <class SCALAR, class PART, class REF>
  void stoch_imag_time_detail(const pcps::Input & userinp,
                              const pcps::Hamiltonian<SCALAR> & ham,
                              pcps::Wavefunction<SCALAR> & wfn);

  //--------------------------------------------------------------------------------
  // pcps::update_correlator -- updates a single correlator by diagonalizing a local
  //                            Hamiltonian
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void update_correlator(const pcps::Input & userinp,
                         const pcps::Hamiltonian<SCALAR> & ham,
                         const pcps::Wavefunction<SCALAR> & wfn,
                         pcps::Correlator<SCALAR> & corr,
                         const double noise_magnitude,
                         const double damping,
                         const int perturbation_level);

  //-------------------------------------------------------------------------------
  // pcps::hartree_fock_spinless -- computes the 1-body reduced density matrix of
  //                                the spinless Hartree-Fock wavefunction
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void hartree_fock_spinless(const pcps::Input & userinp,
                             const pcps::Hamiltonian<SCALAR> & ham,
                             std::vector<SCALAR> & rdm1);

  //-------------------------------------------------------------------------------
  // pcps::hartree_fock_restricted -- computes the 1-body reduced density matrix of
  //                                  the restricted Hartree-Fock wavefunction
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void hartree_fock_restricted(const pcps::Input & userinp,
                               const pcps::Hamiltonian<SCALAR> & ham,
                               std::vector<SCALAR> & restricted_orb_coeffs);

  //-------------------------------------------------------------------------------
  // pcps::hartree_fock_restricted_quad -- computes the 1-body reduced density
  //                                       matrix of the restricted Hartree-Fock
  //                                       wavefunction using the quadratically
  //                                       convergent Newton-Raphson algorithm
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void hartree_fock_restricted_quad(const pcps::Input & userinp,
                                    const pcps::Hamiltonian<SCALAR> & ham,
                                    std::vector<SCALAR> & rdm1);

  //-------------------------------------------------------------------------------
  // pcps::free_particle_ref -- computes the 1-body reduced density matrix of the
  //                            free particle wavefunction
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void free_particle_ref(const pcps::Input & userinp,
                         const pcps::Hamiltonian<SCALAR> & ham,
                         std::vector<SCALAR> & rdm1);

  //--------------------------------------------------------------------------------
  // pcps::correlation_functions -- computes some correlation functions
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void correlation_functions(const pcps::Input & userinp,
                             const pcps::Wavefunction<SCALAR> & wfn);

  template <class SCALAR, class PART, class REF>
  void correlation_functions_detail(const pcps::Input & userinp,
                                    const pcps::Wavefunction<SCALAR> & wfn);

  //-------------------------------------------------------------------------------
  // pcps::type_lookup -- a function that returns the name of the type whose info
  //                      is given
  //-------------------------------------------------------------------------------

  std::string type_lookup(const std::type_info & info);

  //-------------------------------------------------------------------------------
  // pcps::send_vector -- Send a vector to another process
  //                      The size of the vector's data type T must not be variable
  //-------------------------------------------------------------------------------

  template <class T> inline void send_vec(const MPI::Comm & comm,
                                          const std::vector<T> & vec,
                                          const int dest,
                                          const int tag)
  {

    // determine the size of the vector
    const int size = sizeof(T) * vec.size();

    // send the size of the vector
    comm.Send((const void *)(&size), 1, MPI::INT, dest, tag);

    // if the vector is not empty, send its elements
    if (size > 0)
      comm.Send((const void *)(&vec[0]), size, MPI::CHAR, dest, tag);

  }

  //-------------------------------------------------------------------------------
  // pcps::recv_vector -- Receive a vector from another process
  //                      The size of the vector's data type T must not be variable
  //-------------------------------------------------------------------------------

  template <class T> inline void recv_vec(const MPI::Comm & comm,
                                          std::vector<T> & vec,
                                          const int source,
                                          const int tag)
  {

    // receive the size of the vector
    int size;
    comm.Recv((void *)(&size), 1, MPI::INT, source, tag);

    // resize the vector
    vec.resize(size / sizeof(T));

    // if the vector is not empty, receive its elements
    if (size > 0)
      comm.Recv((void *)(&vec[0]), size, MPI::CHAR, source, tag);

  }

  //-------------------------------------------------------------------------------
  // pcps::allocate_shared_array -- allocate an array using a shared pointer
  //-------------------------------------------------------------------------------

  template <class T>
  inline boost::shared_array<T> allocate_shared_array(const int n) {
    
    // check that the length of the array is non-negative
    assert(n >= 0);

    // allocate the array
    boost::shared_array<T> retval(new T[n]);

    // return a shared pointer to the array
    return retval;
  }

  //-------------------------------------------------------------------------------
  // pcps::allocate_shared_array -- allocate and initialize an array using a
  //                                shared pointer
  //-------------------------------------------------------------------------------

  template <class T>
  inline boost::shared_array<T> allocate_shared_array(const int n, const T & value) {

    // allocate the array
    boost::shared_array<T> retval = pcps::allocate_shared_array<T>(n);

    // initialize all the elements to the supplied value
    T * ptr = retval.get();
    T * end_ptr = retval.get() + n;
    while (ptr != end_ptr) {
      *ptr = value;
      ptr++;
    }

    // return a shared pointer to the array
    return retval;
  }

  //-------------------------------------------------------------------------------
  // pcps::unity -- returns the number one in the appropriate type
  //-------------------------------------------------------------------------------

  // no specialization
  template <class SCALAR> inline SCALAR unity() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::unity() for type \"%s\"")
                             % typeid(SCALAR).name() ).str() );
  }
  // specialization for double
  template <> inline double unity<double>() { return 1.0; }
  // specialization for complex
  template <> inline std::complex<double> unity< std::complex<double> >() { return std::complex<double>(1.0, 0.0); }

  //-------------------------------------------------------------------------------
  // pcps::zero -- returns the number zero in the appropriate type
  //-------------------------------------------------------------------------------

  // no specialization
  template <class SCALAR> inline SCALAR zero() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::zero() for type \"%s\"")
                             % typeid(SCALAR).name() ).str() );
  }
  // specialization for double
  template <> inline double zero<double>() { return 0.0; }
  // specialization for complex
  template <> inline std::complex<double> zero< std::complex<double> >() { return std::complex<double>(0.0, 0.0); }

  //-------------------------------------------------------------------------------
  // pcps::imaginary_unity -- returns the number i in the appropriate type
  //-------------------------------------------------------------------------------

  // no specialization
  template <class SCALAR> inline SCALAR imaginary_unity() {
    throw pcps::Exception( ( boost::format("no specialization of pcps::imaginary_unity() for type \"%s\"")
                             % typeid(SCALAR).name() ).str() );
  }
  // specialization for double
  template <> inline double imaginary_unity<double>() { return 0.0; }
  // specialization for complex
  template <> inline std::complex<double> imaginary_unity< std::complex<double> >() { return std::complex<double>(0.0, 1.0); }

  //-------------------------------------------------------------------------------
  // pcps::format_number -- formats a number into a string
  //-------------------------------------------------------------------------------

  template <class SCALAR> inline std::string format_number(const std::string & f, const SCALAR value) {
    return (boost::format(f) % value).str();
  }
  // specialization for complex
  template <> inline std::string format_number< std::complex<double> >(const std::string & f, const std::complex<double> value) {
    std::string retval;
    retval.append("( ");
    retval.append( (boost::format(f) % value.real()).str() );
    retval.append(", ");
    retval.append( (boost::format(f) % value.imag()).str() );
    retval.append(" )");
    return retval;
  }

  //-------------------------------------------------------------------------------
  // pcps::conj -- returns the complex conjugate of a number
  //-------------------------------------------------------------------------------

  template <class SCALAR> inline SCALAR conj(const SCALAR value) {
    return value;
  }
  // specialization for complex
  template <> inline std::complex<double> conj< std::complex<double> >(const std::complex<double> value) {
    return std::conj(value);
  }
  template <> inline std::complex<long double> conj< std::complex<long double> >(const std::complex<long double> value) {
    return std::conj(value);
  }
  template <> inline std::complex<float> conj< std::complex<float> >(const std::complex<float> value) {
    return std::conj(value);
  }

  //-------------------------------------------------------------------------------
  // pcps::real -- returns the real part of a number
  //-------------------------------------------------------------------------------

  template <class SCALAR> inline double real(const SCALAR value) {
    return double(value);
  }
  // specialization for complex
  template <> inline double real< std::complex<double> >(const std::complex<double> value) {
    return std::real(value);
  }

  //-------------------------------------------------------------------------------
  // pcps::imag -- returns the imaginary part of a number
  //-------------------------------------------------------------------------------

  template <class SCALAR> inline double imag(const SCALAR value) {
    return 0.0;
  }
  // specialization for complex
  template <> inline double imag< std::complex<double> >(const std::complex<double> value) {
    return std::imag(value);
  }

  //-------------------------------------------------------------------------------
  // pcps::square_norm -- returns the square of the norm of a number
  //-------------------------------------------------------------------------------

  inline double square_norm(const double x) {
    return x*x;
  }
  inline double square_norm(const std::complex<double> x) {
    return x.real() * x.real() + x.imag() * x.imag();
  }

  //-------------------------------------------------------------------------------
  // pcps::round -- rounds a double to the nearest integer
  //-------------------------------------------------------------------------------

  inline int round(const double d) {
    return ( d >= 0 ? int(d+0.5) : int(d-0.5) );
  }

  //-------------------------------------------------------------------------------
  // pcps::round -- rounds a number to a given digit precision
  //-------------------------------------------------------------------------------

  inline double round(const double x, const int digits) {
    const double factor = std::pow(double(10.0), digits);
    const double y = factor * x;
    return double( y >= 0.0 ? (long int)(y+0.5) : (long int)(y-0.5) ) / factor;
  }

  inline std::complex<double> round(const std::complex<double> x, const int digits) {
    return std::complex<double>( pcps::round(x.real(), digits), pcps::round(x.imag(), digits) );
  }

  //-------------------------------------------------------------------------------
  // pcps::is_half_integer -- checks whether a double is equal to some integer
  //                          divided by two
  //-------------------------------------------------------------------------------

  inline bool is_half_integer(const double x) {
    return ( std::fabs( pcps::round(2.0 * x) - 2.0 * x ) < 1.0e-9 );
  }

  //-------------------------------------------------------------------------------
  // pcps::get_classical_tripartite_guess -- construct an initial guess for a
  //                                         tripartite lattice in which each
  //                                         sublattice's spin's expectation value
  //                                         is rotated at an angle of 2*pi/3
  //                                         (120 degrees) from the others
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void get_classical_tripartite_guess(const pcps::Input &userinp,
                                      const std::set<int> & a_site_set,
                                      const std::set<int> & b_site_set,
                                      const std::set<int> & c_site_set,
                                      std::vector< pcps::Correlator<SCALAR> > &correlators);

  //-------------------------------------------------------------------------------
  // pcps::get_classical_icosidodecahedron_guess -- construct an initial guess
  //                                                using single site correlators
  //                                                to express the spin coherent
  //                                                state representation of the
  //                                                classical heisenberg solution
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void get_classical_icosidodecahedron_guess(const pcps::Input &userinp, std::vector< pcps::Correlator<SCALAR> > &correlators);

  //-------------------------------------------------------------------------------
  // pcps::get_classical_triangular_guess -- construct an initial guess using
  //                                         single site correlators to express the
  //                                         spin coherent state representation of
  //                                         the classical heisenberg solution of
  //                                         the triangular lattice
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void get_classical_triangular_guess(const pcps::Input &userinp, std::vector< pcps::Correlator<SCALAR> > &correlators);

  //-------------------------------------------------------------------------------
  // pcps::nr_get_residual -- computes the newton raphson residual vector
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void nr_get_residual(const pcps::Input & userinp,
                       const pcps::Hamiltonian<SCALAR> & ham,
                       const pcps::Wavefunction<SCALAR> & wfn,
                       std::vector<SCALAR> & residual);

  //-------------------------------------------------------------------------------
  // pcps::nr_get_jacobian -- computes the newton raphson jacobian matrix
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void nr_get_jacobian(const pcps::Input & userinp,
                       const pcps::Hamiltonian<SCALAR> & ham,
                       const pcps::Wavefunction<SCALAR> & wfn,
                       std::vector<SCALAR> & jacobian);

  //-------------------------------------------------------------------------------
  // pcps::projected_nr_update -- updates the correlators via a Newton-Raphson step
  //-------------------------------------------------------------------------------

  template <class SCALAR>
  void projected_nr_update(const pcps::Input & userinp,
                           const pcps::Hamiltonian<SCALAR> & ham,
                           pcps::Wavefunction<SCALAR> & wfn,
                           const int perturbation_level,
                           const double damping);

  //--------------------------------------------------------------------------------
  // pcps::solve_projected_nr_equations -- Solves the newton raphson equations
  //                                       of projected CPS for the update vector,
  //                                       which is returned in the residual.
  //                                       The jacobian data may be destroyed.
  //--------------------------------------------------------------------------------

  template <class SCALAR>
  void solve_projected_nr_equations(const pcps::Input & userinp,
                                    const int n,
                                    SCALAR * jacobian,
                                    SCALAR * residual);

  //-------------------------------------------------------------------------------
  // pcps::bcast -- broadcasts data from one process to others
  //-------------------------------------------------------------------------------

  template<class T>
  inline void bcast(const MPI::Comm & comm,
                    T * data,
                    size_t n,
                    const MPI::Datatype & datatype,
                    int root) {
    const int block_size = 100000000;
    while ( n > size_t(block_size) ) {
      comm.Bcast((void *)data, block_size, datatype, root);
      data += block_size;
      n -= size_t(block_size);
    }
    if (n > 0)
      comm.Bcast((void *)data, int(n), datatype, root);
  }

  //-------------------------------------------------------------------------------
  // pcps::reduce -- reduce data to root process
  //-------------------------------------------------------------------------------

  inline void reduce(const MPI::Comm & comm, const double * send_buf, double * recv_buf, int n, const MPI::Op & op, int root) {
      comm.Reduce(send_buf, recv_buf, n, MPI::DOUBLE, op, root);
  }
  inline void reduce(const MPI::Comm & comm, const std::complex<double> * send_buf, std::complex<double> * recv_buf, int n, const MPI::Op & op, int root) {
      comm.Reduce(send_buf, recv_buf, n, MPI::DOUBLE_COMPLEX, op, root);
  }
  inline void reduce(const MPI::Comm & comm, const int * send_buf, int * recv_buf, int n, const MPI::Op & op, int root) {
      comm.Reduce(send_buf, recv_buf, n, MPI::INT, op, root);
  }
  inline void reduce(const MPI::Comm & comm, const long int * send_buf, long int * recv_buf, int n, const MPI::Op & op, int root) {
      comm.Reduce(send_buf, recv_buf, n, MPI::LONG, op, root);
  }

  //-------------------------------------------------------------------------------
  // pcps::allreduce -- reduce data to all processes
  //-------------------------------------------------------------------------------

  inline void allreduce(const MPI::Comm & comm, const double * send_buf, double * recv_buf, int n, const MPI::Op & op) {
      comm.Allreduce(send_buf, recv_buf, n, MPI::DOUBLE, op);
  }
  inline void allreduce(const MPI::Comm & comm, const std::complex<double> * send_buf, std::complex<double> * recv_buf, int n, const MPI::Op & op) {
      comm.Allreduce(send_buf, recv_buf, n, MPI::DOUBLE_COMPLEX, op);
  }
  inline void allreduce(const MPI::Comm & comm, const int * send_buf, int * recv_buf, int n, const MPI::Op & op) {
      comm.Allreduce(send_buf, recv_buf, n, MPI::INT, op);
  }
  inline void allreduce(const MPI::Comm & comm, const long int * send_buf, long int * recv_buf, int n, const MPI::Op & op) {
      comm.Allreduce(send_buf, recv_buf, n, MPI::LONG, op);
  }

  //-------------------------------------------------------------------------------
  // pcps::strip_comments -- removes all text after a comment character
  //-------------------------------------------------------------------------------

  inline void strip_comments(std::string & s, const char comment_char) {
    for (std::string::iterator c = s.begin(); c != s.end(); c++)
      if (*c == comment_char) {
        s.erase(c, s.end());
        break;
      }
  }

  //-------------------------------------------------------------------------------
  // pcps::cmpd_pair_index -- Computes a compound index for a pair of distinct
  //                          numbers i < j in the range 0, 1, 2, ..., n-1.
  //-------------------------------------------------------------------------------

  inline int cmpd_pair_index(const int i, const int j, const int n) {
    assert( i < j );
    assert( i >= 0 && i < n );
    assert( j >= 0 && j < n );
    return ( 2*n - i - 3 ) * i / 2 + j - 1;
  }

}

#endif
